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Abstract 

We present results from numerical simulations of three different 3d four-fermion 
models that exhibit Z 2 , U(l), and 577(2) x 577(2) chiral symmetries, respec- 
tively. We performed the simulations by using the hybrid Monte Carlo algorithm. 
We employed finite size scaling methods on lattices ranging from 8 3 to 40 3 to 
study the properties of the second order chiral phase transition in each model. 
The corresponding critical coupling defines an ultraviolet fixed point of the renor- 
malization group. In our high precision simulations, we detected next-to-leading 
order corrections for various critical exponents and we found them to be in good 
agreement with existing analytical large- Nf calculations. 
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1 Introduction 



The 3d four-fermion models are among the simplest relativistic quantum field 
theories of interacting fermions. There are several motivations for studying such 
models. Dynamical breaking of chiral symmetry occurs at strong enough interac- 
tion coupling g 2 c . The chirally broken phase is separated from the chirally sym- 
metric phase by a second order phase transition at the critical coupling. Even 
though these models are not perturbatively renormalizable, it has been shown 
that the 1/Nf expansion about the fixed point g\ is exactly renormalizable [Q]|. 
In addition, four-fermion models are ideal laboratories for studying continuum 
phase transitions in the presence of massless fermions. Hence, they define new 
universality classes that are quantitatively different from the ferromagnetic phase 
transitions in bosonic O(N) Heisenberg spin models. Furthermore, in the frame- 
work of 1/Nf expansion 0, it has been shown that the universality class of the 
d-dimensional four-fermion models, where d is between two and four, is the same 
as the universality class of the Higgs-Yukawa model with the same chiral sym- 
metry. Understanding the properties of the continuum phase transition, which 
separates the chirally symmetric from the chirally broken phase, requires non- 
perturbative techniques such as the large- Nf expansion Iffll3ll4l [51161171 l8ll9l[T0ll. 
exact renormalization group equations [fTTl[T2l . and lattice Monte Carlo simula- 
tions El [HO. 

Given that 3d four-fermion models incorporate certain important features of 
QCD, they have been used recently as model field theories to study the proper- 
ties of the strong interaction at non-zero temperature and non-zero quark number 
density [15). In addition, there may be applications of four-fermion models to 
high-T c superconductivity |fT6ll , for instance in describing non-Fermi liquid be- 
havior in the normal phase ifPTTl . More recently, it was proposed that the Hubbard 
model on a honeycomb lattice relevant to the newly discovered graphene sheets, 
has a transition described by the three-dimensional Z 2 -symmetric four-fermion 
model fig! . 

In this paper, we study numerically the critical properties of the 3c? four- 
fermion models that exhibit the three different Z 2 , abelian 17(1), and non-abelian 
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SU (2) x SU (2) chiral symmetries. In our simulations, we fixed the number of 
fermion flavors to Nf = 4. Our simulations are the first accurate finite size scaling 
(FSS) studies of the 17(1) and SU{2) x SU(2) models that allow us to detect next- 
to-leading order corrections on the values of the critical exponents and to compare 
them with existing analytical large- Nf predictions [|51 |6j Ul [HI 13 QZ3- Our results 
from the Z 2 model simulations are also in good agreement with existing large-iVj 
predictions as are other accurate Monte Carlo results with Nf = 2 lfi~3l . 



2 Models and Observables 

In this section, we introduce the three different versions of the model we shall 
be dealing with in the bulk of this paper and the observables used to measure 
the critical exponents of the continuous phase transitions. In the literature, the 
models are often called the Gross-Neveu models and their continuum space-time 
lagrangians (we work in Euclidean space) are as follows: 

L A = +m)*i- i^-i^i^i? (1) 
L B = ViW + m)^ - 2^-[(*i*i) 2 - (*i75*<) 2 ]- (2) 

C c = f i(0 +m)Vi-£ [(Vi%) 2 - {%jst%) 2 ] . (3) 

We treat as four-component Dirac spinors and the index i runs over Nf 

fermion species. It can be easily shown that in the chiral limit m — > 0, La has a 
Z 2 , L B a U(l), and L c an SU{2) x SU{2) chiral symmetry. 

For analytical and computational purposes, it is useful to introduce auxiliary 
fields a and 7Tj. Hence, the bosonized lagrangians become quadratic in tyf. 



Nt 

^ = ^(^+m + a)^ + -4(T 2 . (4) 

^9 



Nf 

L B = ^ l (^+m + a + i l5 7r)^ l + ^(a 2 + n 2 ) (5) 
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C c = #i + m + a + Z757T ■ f) % + — | (a 2 + t? ■ tt) . (6) 

For sufficiently strong coupling g 2 > g 2 the models exhibit spontaneous symmetry 
breaking implying dynamical generation of a fermion mass. The pion fields 7Tj 
become the associated Goldstone bosons. 

We used the staggered fermion discretization with the auxiliary fields living on 
the dual lattice sites to formulate the models in their bosonized form on the lattice. 
For each case, we used the hybrid Monte Carlo algorithm with Nf = 4 fermion 
flavors to perform numerical simulations exactly. The Monte Carlo procedure 
was optimized by choosing the microcanonical trajectory length at random from 
a Poisson distribution with mean value equal to 1.0. This method of optimization 
which guarantees ergodicity was found to decrease the autocorrelations in the 
data significantly |fi~9l . Details concerning the lattice actions and the numerical 
algorithm can be found in Il4~l l20"ll211 . 

We work in the chiral limit to study the chiral phase transition of the models. 
Hence, we choose not to introduce a bare quark mass into the lattice action. With- 
out the benefit of this interaction, the direction of symmetry breaking changes over 
the course of the simulation such that E = y J2 X a ( x ) an( ^ ^k = y J2 X n i( x ) av ~ 
erage to zero over the ensemble. It is in this way that the absence of spontaneous 
symmetry breaking on a finite lattice is enforced. Another option is to introduce 
an effective order parameter $ equal to the magnitude of the vector $ = (E, II). 
In the thermodynamic limit, ($) is equal to the true order parameter (a) extrapo- 
lated to zero quark mass. 

We employ the finite size scaling (FSS) method [221, a well-established tool, 
to study the critical behavior of the model on lattices available to us. The correla- 
tion length £ on a finite lattice is limited by the size of the system and consequently 
no true criticality can be observed. The dependence of a given thermodynamic ob- 
servable, A, on the size L of the box is singular. According to the FSS hypothesis, 
in the large volume limit, A is given by: 

A(t,L)=L^Q A (tL 1 /»), (7) 
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where t = {(3 C — (3)/ [3 C is the reduced temperature, u is the exponent of the cor- 
relation length, Qa is a scaling function that is not singular at zero argument, and 
Pa is the critical exponent for the quantity A. Using eq. ©, one can determine 
such exponents by measuring A for different values of L. 

In the large L limit, the FSS scaling form of the effective order parameter ($) 
is given by 

($) = L~ pm/u f a (tL l/u ). (8) 

A standard method to measure the inverse critical coupling [3 C = 1/g 2 for a 
second order transition is to compute the Binder cumulant U B ((3, L) [|23l , defined 
by 

1 ($ 4 ) 

for various system sizes. Near the critical coupling and on sufficiently large lat- 
tices, where subleading corrections from the finite lattice size L are negligible, 
Ub = fBL(tL 1 ' v ). Therefore, at (3 C , Ub becomes independent of L. Deviations 
from this relation can be explained by finite size confluent corrections. The lead- 
ing Li/L dependence in the deviation of the intersection point j3* from the critical 
point [3 C is estimated by Binder ll23l as 
1 1 a 

~^{L)~J c + \n{L 1 /Ly ( } 

In our analysis we chose L to be the smallest lattice size L = 8 and hence L\ are 
the remaining lattice sizes. 

For the general 0(n) -symmetric models, it can be easily shown ll24ll that as the 
lattice volume tends to infinity in the weak coupling limit, Gaussian fluctuations 
around $ = lead to Ub — ► 2(n — l)/3n. For n — 1 (0(1) = Z 2 symmetry) this 
gives a zero reference point, for n = 2 (0(2) = U(1) symmetry) U B — » 1/3, and 
for n = 4 (0(4) = SU(2) x SU(2) symmetry) U B -> 1 /2. In the chirally broken 
phase Ub — > 2/3 for all n in the thermodynamic limit. 

Another quantity of interest is the susceptibility x that is given, in the static 
limit of the fluctuation-dissipation theorem, by 

X = lim V[(&) - ($) ■ ($)], (11) 
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where V is the lattice volume. For finite systems, the true order parameter (<$>) 
vanishes and for (3 > (3 C the susceptibility is given by: 

X = V{&). (12) 

This observable should scale at criticality like 

X = L^f x (tL^). (13) 

Furthermore, the logarithmic derivatives of ($) can give estimates for the crit- 
ical exponent v. It can be easily shown that 



(S b ) 



(14) 



(<&} 

where S b is the bosonic part of the lattice action that is multiplied by the coupling 
(3. D has a scaling relation 

D = L^f D (tL^)- (15) 

We used the histogram reweighting method E51 to perform our study most 
effectively. This enabled us to calculate the observables in a region of couplings 
around the simulation coupling. We utilized this technique efficiently by perform- 
ing simulations at slightly different couplings $ close to the critical coupling (3 C . 
We also employed the jacknife method to estimate the statistical errors on the var- 
ious observables reliably. This method accounts for correlations in the raw data 
set. 



3 Results 

In this section we present the results of the data analysis for the three different 
models. In all three cases, the fermion species number is fixed at Nf = 4. An 
accurate determination of the critical exponents requires a precise determination 
of the critical coupling. We calculated the critical couplings by using the Binder 
cumulant technique described in the previous section. For different lattice sizes, 
the curves U B = U B {(3) should intersect at (3 = (3 C up to finite size corrections that 
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Figure 1: Binder cumulant vs. (3 for different lattice sizes; Z 2 model. 

are visible on the smaller lattices. We used the histogram reweighting method to 
obtain the values of Ub versus (3. We show these values for the Z 2 model in Fig.Q] 
We performed the simulations on the largest 40 3 lattice at a single value of the 
coupling j3 = 0.835 and we generated approximately half a million configurations 
with average trajectory length equal to 1 .0. We performed the simulations on the 
other lattices (8 3 , 12 3 , 16 3 , 22 3 , 30 3 ) at all values (3 = 0.82, 0.83, 0.84, 0.85 with 
approximately half to one million configurations for each (3. It is clear that in the 
Z 2 model, the U B curves intersect at (f3 c , U B {f3 c )) = (0.835(1), 0.232(8)). 

As expected, the situation is somewhat different in the U(l) model. In this 
model infrared fluctuations are stronger than in the discrete symmetry model. As 
a result, finite size effects near the critical coupling are larger for U(l) than for 
Z 2 . We performed the simulations for the U(l) model on the 40 3 lattices at all 
values (3 = 0.830, 0.835, 0.840, 0.845, 0.850, 0.86, whereas on the smaller lattices 
at all values f3 = 0.83, 0.86 and in steps of 0.01. The data set generated on 
30 3 at P = 0.850 was corrupted and it was not included in the analysis. Approx- 
imately 6 x 10 5 - 1.3 x 10 6 configurations were generated at each (3. We show 
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the values of Ub versus (3 in Fig. [2] The leading L\/L finite size corrections are 
taken into consideration by using eq. ( TTOT ). We plot (l/ln(Li/L), for L = 8 
in Fig. [3] We computed the errors for from the jacknife errors for Ub(P)- 
The extrapolation of jjj^ to the point 1/ m(Li/L) = gives (3 C = 0.853(2) and 
Ub(Pc) = 0.424(8) on the 40 3 lattice. 

We performed an analysis for the U b data for SU (2) x SU(2) similar to the one 
forthe£/(l) model. In this case, we performed simulations at (3 = 0.92, 1.00 in 
steps of 0.02 for the 8 3 , 12 3 , 16 3 , 22 3 , 30 3 lattices and at p = 0.94, 0.98 in steps 
of 0.01 on the largest 40 3 lattice. The curves U b versus j3 obtained from histogram 
reweighting at two consecutive values of j3 did not intersect. Therefore, to obtain 
the intersection we used a linear approximation in the middle region between two 
curves. The values of Ub on different lattices near (3 C are shown in Fig. |4]and 
the extrapolation of 1//?* to the point \/hx(L\/ L) = are shown in Fig.|3j We 
extracted from this analysis the values (3 C = 0.960(3) and U(j3 c ) = 0.544(7) on 
the largest 40 3 lattice. 




8 



1.3 

1.25 



/3* 



k/(2) x sdm 




HLi/L) 



Figure 3: The intersection of Ub(L) and Ub{L\) for L = 8 vs. ln(Li/L). 
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Figure 5: Effective order parameter ($) as a function of the lattice size L for all 
three models. 




Next, we calculated the exponent ratios ft m /v for the three models by fit- 
ting to eq. © the values of (<E>) at (3 C obtained on different lattice sizes. After 
fitting the data obtained on all lattice sizes we get (3 m /u = 0.927(15) for Z 2 , 
Pm/v = 0.955(20) for 17(1), and fi m /v = 1.04(2) for SU(2) x SU(2). These 
values of (3 m /u take into consideration the statistical error in [3 C . In an effort to 
check to what extend our results are affected by possible small volume effects we 
repeated the analysis without including the smallest lattice. Our results, summa- 
rized in table Q] show that any finite volume systematic effects are smaller than the 
statistical errors. Another analysis where the 12 3 data were excluded confirmed 
this conclusion. The data and the fitted functions (for L = 12, ...40) for the three 
models are shown in Fig. [5] 

Similarly, we obtained the exponent ratios j/v by fitting the data for the sus- 
ceptibility x (eq- (©) at f3 c to its FSS relation eq. (flSb . We present our results in 
Table [2] together with analytical predictions obtained from large- Nj calcultations 
to order 1 /Nj [H. It is clear that our numerical results are in good agreement 
with the analytical predictions. Furthermore, the results we got after omitting the 
smallest 8 3 volume show that any finite size systematic effects are within the sta- 
tistical errors. The data and the fitted functions (for L = 12, 40) for the three 
models are shown in Fig. [6] 

We used the logarithmic derivative D, defined in eq. (fl~4~l) . to calculate the 
exponent v. According to eq. (fT5l) . at f3 c , D ~ L 1 ^ . We present the values of v 
for each model in Table [3] together with the respective values obtained from large- 
Nf calculations to order 1/Nj ||71|9l[T0|. As in the j/u case, the results obtained 
from our simulations are in good agreement with the analytical predictions. Like 
in the previous observables, in this case also systematic small volume effects are 
within the statistical errors. The data and the fitting functions (for L = 12, ...,40) 
for the three models at their critical couplings are shown in Fig. [7] 

Using our results for /3 m /u and 7///, obtained from fits on all lattice sizes, we 
can check whether the hyperscaling relation 

^ + 12-^ = (16) 
v 1v 2 

is satisfied. We find that for all three different models the left hand side of eq. (fTBI) 
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Table 1: Values of j3 m /v measured from our simulations. 





z 2 


U(l) 


SU{2) x SU(2) 


simulations L = 8, 40 


0.927(15) 


0.955(20) 


1.04(2) 


simulations L = 12, 40 


0.917(20) 


0.952(25) 


1.05(3) 



is consistent with the value zero with a statistical uncertainty of 3-4%. 




Figure 7: Logarithmic derivarive D of the order parameter as a function of the 
lattice size L for all three models. 



4 Conclusions 

We presented results from Monte Carlo simulations of 3d four-fermion models 
with Z 2 , U(l), and SU(2) x 577(2) chiral symmetries. These models are among 
the simplest relativistic field theories of interacting fermions, and therefore are 
benchmarks for studying critical phenomena in the presence of massless fermions. 
They are also used as model field theories to study the behavior of strong interac- 
tion under extreme conditions and have applications in condensed matter systems. 
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Table 2: Values ofj/u measured from our simulations and from large- Nf calcu- 
lations. 







z 2 


U{\) 


SU{2) x SU{2) 


simulations L = 8, .. 


.,40 


1.152(25) 


1.09(3) 


0.925(25) 


simulations L = 12, . 


..,40 


1.165(40) 


1.09(4) 


0.910(30) 


large iV, 




1.132 


1.06 


0.946 



Table 3: Values v measured from our simulations and from large- Nf calculations. 





z 2 


17(1) 


SU{2) x SU{2) 


simulations L = 8, 40 


0.98(2) 


1.05(2) 


1.14(3) 


simulations L = 12, 40 


0.99(2) 


1.03(4) 


1.16(5) 


large N f H1B\M 


0.98 


1.02 


1.11 



In all three cases, we performed simulations with Nf = 4 fermion species. Ana- 
lytical calculations predict small next-to-leading order corrections for the critical 
exponents of the second order phase transitions of these models at this intermedi- 
ate value of Nf. We detected these corrections in our simulations by employing 
standard finite size scaling techniques and we found them to be in good agreement 
with large- Nf expansions up to 0(1/Nj) j6l|Vl[8l|9l[T0l. Our results improve sig- 
nificantly previous numerical studies of 3d four-fermion models. Future work 
with much better statistics on a variety of lattices including larger sizes than the 
ones used in this work will allow for the detection of corrections to scaling effects 
and possible deviations from the 0(1/ Nj) analytical calculations. Also simula- 
tions with Nf = 1 will be particularly instructive, as for such a small Nf large- Nf 
calculations cannot be applied. 
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